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Abstract 

Using a two-dimensional fluid description, we investigate the nonlinear radial-longitudinal dy- 
namics of intense beams in storage rings and cyclotrons. With a multiscale analysis separating the 
time scale associated with the betatron motion and the slower time scale associated with space- 
charge effects, we show that the longitudinal-radial vortex motion can be understood in the frame 
moving with the charged beam as the nonlinear advection of the beam by the E x B velocity field, 
where E is the electric field due to the space charge and B is the external magnetic field. This 
interpretation provides simple explanations for the stability of round beams and for the develop- 
ment of spiral halos in elongated beams. By numerically solving the nonlinear advection equation 
for the beam density, we find that it is also in quantitative agreement with results obtained in PIC 
simulations. 



ccrfon@cims.nyu.edu 
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I. INTRODUCTION 



In recent years, novel uses of particle accelerators for nuclear energy, nuclear security, 
medicine and material science applications have triggered the design, development and con- 
struction of moderate size, high intensity cyclotrons. The higher beam intensities expected 
in these new machines provide more powerful tools for these new endeavours, but also come 
with very stringent uncontrolled beam loss requirements. Satisfying these low beam loss 
criteria requires detailed knowledge of the beam dynamics in such machines. In particular, 
collective modes and instabilities associated with space charge effects have to be taken into 
account and need to be studied and understood. Theory and simulation tools are being 
developed to that purpose. The computational effort is crucial because of the complexity 
of the problem at hand. One can hardly imagine a purely analytic study that would accu- 
rately account for the complex magnetic geometry of modern cyclotrons and at the same 
time properly model the nonlinear beam dynamics due to space charge effects. However, 
analytical theory work is equally as important in order to provide an interpretation for the 
results from complicated, all-inclusive simulations as well as to identify the basic contributing 
mechanisms. 

Interestingly, the computational effort in high-intensity cyclotrons is almost entirely fo- 
cused on Particle-In-Cell (PIC) methods [lj] , 3| , , Q] . PIC methods have a number of ad- 
vantages. They are intuitive, naturally lead to parallel solvers, and can easily take exact 
magnetic geometries as inputs. Nowadays, PIC methods are commonly used to design, in- 
terpret and corroborate experiments in which space charge effects play a significant role in 
the beam dynamics Q|,0|,j7]. However, precisely because of their conceptual simplicity, PIC 
codes tend not to yield as much insights into the physical mechanisms at work in a given 
observed phenomena as do continuum descriptions. Theoretical work based on continuum 
descriptions, whether kinetic or fluid, are often better suited to identify the key nondimen- 
sional parameters and to derive the typical scalings for the dependence of the quantities 
of interest on these nondimensional parameters. Analytic and numerical work based on 
continuum descriptions are therefore a necessary and useful complement to PIC simulations 

In this article, we illustrate the constructive interplay of experiments, PIC simulations and 
theory in the study of a topic of high interest in beam dynamics in isochronous cyclotrons, 
namely the coherent longitudinal-radial vortex motion associated with the nonlinear inter- 
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olayof radial and longitudinal space-charge forces. It has been oberved in PIC simulations 
3|;|8|, 9], 10( that a coasting beam which is elongated in the longitudinal direction (the direc- 
tion of propagation of the beam) develops a spiral galaxy shape due to space-charge effects. 
This phenomenon, which has sometimes been called the "spiraling instability" can play 



a positive role when the spiral arms are lost after a few turns, at 



center of the spiral becomes a very compact stable beam 



3,0 



ow beam energy, and the 



121 ] . The spiraling of 



the beam can also play a negative role. It is thought to lead to beam breakup, as confirmed 
by experiment and PIC simulations [l3j],[l4j]. At lower beam currents, the spiral arms can 
also grow over longer time scales, and therefore form a beam halo only after a large number 
of turns, at which point it has been accelerated to high energies. 

Until now, aside from PIC simulations, the theoretical explanations of beam spiraling 
have bee. based oa single-particle pictures Q.fiQ.Q. ba these single-particle de- 
scriptions, the evolution of the beam in the presence of space-charge forces is extrapolated 
from the modifications of the motion of single particles due to this force. In contrast, using 
a fluid model of the charged beam we present a new self-consistent nonlinear description 
of the coherent radial-longitudinal vortex motion and give an intuitive explanation for the 
formation of the spiral galaxy shape in an elongated charged particle beam. With a mul- 
tiscale analysis that separates the time scale associated with the betatron motion and the 
time scale associated with the space-charge force, we show that in the frame rotating with 
the center of the beam, the evolution of the beam on the space-charge time scale is deter- 
mined by the nonlinear advection of the beam density by the E x B velocity. Here, E is 
the self electric field, and B is the external cyclotron magnetic field. With this new fluid 
picture of the radial-longitudinal vortex motion, the observed behavior in experiments and 
PIC simulations is easily understood, and beam spiraling has a simple explanation. 

Consider first a round beam, whose charge density is only a function of its radius. In the 
frame rotating with the beam, the self electric field of such a charge distribution is only along 
the direction of the beam radius, and the E x B velocity vector in the radial-longitudinal 
plane is therefore purely perpendicular to the beam radius at every point. This implies that 
the advection in the E x B field will have no effect on the radial distribution of the beam 
density: a round beam thus has a stationary distribution. In an elongated beam however, 
the cylindrical symmetry of the problem is broken since the density distribution depends 
both on radius and angle. The beam density distribution will thus be distorted by the 
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E x B velocity field. Since the E x B velocity is largest at the extremities of the beam, 
the distortion will be largest away from the beam center, which leads to the formation of 
spiral galactic arms. This intuitive explanation is confirmed by simulations in which we 
numerically solve the nonlinear fluid equations, and we demonstrate quantitative agreement 
between our fluid description and PIC simulations. 

The structure of the paper is as follows. In Section [Til we describe the geometry of the 
problem and introduce the fluid description of the beam which we use for our analysis. 
In Section IIII| we carry out a multiscale expansion of the fluid equations, and derive the 
advection equation for the evolution of the beam density on the time scale associated with 
the space charge force. We solve this equation numerically in Section [TV] for different cases 
of interest, and demonstrate both the formation of spiral galactic arms for elongated beams 
and the stability of round beams. We end with conclusions in Section V. 



II. FLUID DESCRIPTION OF THE BEAM 



One of the main motivations for this work is to offer an analytic derivation of the self- 
consistent beam vortex motion under the effect of space charge. Such a goal implies that we 
will have to make simplifications in order to make our analysis tractable. In making these 
simplifications, our guiding principle is to simplify the geometric aspects of the problem 
enough so that only a minimal amount of beam physics assumptions need to be made. 



A. Geometry of the problem 

It has been observed in PIC simulations that the details of the beam description along 
the vertical direction, i.e. along the direction of the main applied ma gne tic field, have a 
very limited influence on the longitudinal-radial vortex motion Isj],!^!, 16]. To make our 
analysis more tractable, we therefore restrict our study to the two-dimensional problem of a 
nonrelativistic coasting beam in the radial-longitudinal plane, immersed in a homogeneous 
magnetic field along the z axis, which coincides with the vertical direction. Throughout 
our analysis we will thus have B = Be z with B = constant, and d/dz = will hold. One 
implication of these assumptions is that the formation of spiral arms does not depend on 
relativistic effects or field gradient along the vertical direction. 
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In the radial-longitudinal plane, i.e. the plane perpendicular to the magnetic field, it is 
most convenient to study space-charge effects in the frame moving with the coasting beam. 
The transformation from the laboratory frame to the moving frame is best understood in 
terms of the polar coordinates (R, $) associated with the cyclotron geometry and shown in 
Figure [TJ 




FIG. 1: Polar coordinates (R, $) used to define the change of reference frame to the frame rotating 
with the beam, and (in red) Cartesian coordinate system used in the rotating frame 

We consider a single-species intense ion beam, and for the simplicity of the notation, 
the ions are singly charged: q = e. The ion mass is m, so that the cyclotron frequency is 
u c = eB/m. The change of reference frame to the frame rotating with the beam is given by 
the transformation 

R' = R $' = $ + oj c t t' = t v' R = v R v'q = v $ + Rco c (1) 
where the primed quantities represent quantities in the moving frame. Using the relations 
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in Eq.flT]), we easily find 

d d d d d d d 

dR = dR> d¥ = dl 7 di = dt i+ Uc d¥ ^ ' 

In the next section, we will use these relations to express the fluid equations in the frame 
moving with the beam. Before doing so, we briefly discuss the two coordinate systems we 
will use in the local frame moving with the beam, which are shown in Figure [2j The first 
system is a two-dimensional Cartesian coordinate system (x, y) with the origin at the beam 
center, and defined such that (x, y, z) is a right-handed orthogonal coordinate system with 
the x-axis coinciding with the radial direction, i.e the main radius of the cyclotron, and the 
y-axis coinciding with the longitudinal direction. This is the coordinate system also plotted 
in red in Figure [TJ In this coordinate system, the ion beam travels in the — y direction. The 
second system is a polar coordinate system (r, 9) centered at the center of the beam, defined 
such that (r, 9, z) is a right-handed orthogonal cylindrical coordinate system and such that 
x = rcos9 and y = rsin^. It is equipped with the orthonormal basis (e r , eg), where the unit 
vectors e r and are defined in terms of the Cartesian unit vectors and e y according to 
e r = cos^e^. + sin^e^ and eg = — sin^e^. + cos9e y . The Cartesian coordinate system will be 
mostly used for the numerical simulations in section IIVt and the polar coordinate system is 
the most appropriate for the physical interpretation of results. 



B. Fluid description of the beam 

At the beam temperatures and densities relevant to high intensity cyclotrons, the charged 
particle beam is a nonneutral plasma that can be considered collisionless sj]. Then, by taking 
the first two moments of the collisionless Vlasov equation, we obtain the following exact fluid 
equations involving the beam density n, velocity v and pressure tensor P 

^ + V-(nv) = (3) 
fdv \ 

mn I — + v • Vv J + V • P = en (E + v x B) (4) 

In Eq.fllJ), E is the self electric field, which is determined by solving 

en , . 

V ■ E = — 5 
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Radial direction 



FIG. 2: (x, y) and (r, 9) coordinate systems used for the numerical simulations and the physical 
interpretation of results. 

B is the external magnetic field B = Be z since the self magnetic field is negligible in the 



18| . Using the relations in Eq. (jSJ), we can rewrite Equations 



nonrelativistic limit 
and ([5]) in the frame rotating with the beam. We find 



dn' 
~dt' 



+ V ■ (n'v A 







mn' + v' • VV ) + V ■ P' = en' (Ef - v' x B] 



V' -E' 



en 



(6) 
(7) 

(8) 



where as before the primed quantities represent quantities in the moving frame. From now 
on we will always work in the frame rotating with the beam, so to simplify the notation we 
will drop the prime symbols in the remainder of the article. Note that in the moving frame, 
there is a minus in the magnetic field force. It is due to the Coriolis force. 



In the frame rotating with the beam the electrostatic approximation holds {3], so we 
write E = — V</> where 4> is the electrostatic potential. Then, if we call a the characteristic 
size of the beam (its radius if the beam is circular), N its peak density and T its peak 
temperature, renormalizing the fluid quantities according to 



t — > UJ f t 



n — > 



n 



V ->■ aV 
P 



v ->■ 



au r . 



No N T ^ eN a 2T 

leads to the following set of nondimensional fluid equations in the frame rotating with the 



beam: 



dn 

~dt 
dv 

~dt 



+ nV • v = 



v x e. 



-n 



5 2 V0 



a 



VP 



n 



(9) 
(10) 
(11) 



Here, d/dt = d/dt + w- V and 5 2 = uj 2 /u 2 where Co> p is the plasma frequency, uj 2 = N e 2 /meo. 
The parameter a is defined by a 2 = TQ/ma 2 uj 2 = \ 2 D /a 2 where is the Debye length. The 
experimental regime of interest for high intensity beams is a < 1, so that space-charge effects 
dominate over temperature effects. 

It is clear from Eqs.(l9l)- (jTT!) that the parameter determining the importance of space- 
charge effects on the behavior of the beam is 5 2 , the ratio of the beam plasma frequency 
squared to the cyclotron frequency squared. This was observed in a slightly different context 
in Eq.(12) in 13] and in 19j. In particular, in 19j the key role of the parameter x is 
emphasized, with x defined so that x — 5 2 /2 in our notation. As noted in 19j, all cyclotrons 
and rings satisfy 5 2 < 1, and most machines satisfy 5 2 1. The parameter S 2 can then 
be used to derive an approximate form for the pressure tensor P in the limit of small 5, as 
we will do in the next paragraph. It can also be used to perform a multiscale analysis of 
the resulting beam fluid equations which separates the fast motion associated with betatron 
oscillations from the slower motion due to space-charge and thermal forces. This is precisely 
what we do in Section 11111 

Equations (I9l)- (1TT1) cannot be solved as such since we are still facing the problem of 
closure of the hierarchy of moment equations 1^]. However, a key point of our analysis 
is that with our simplified geometry, a single assumption we make concerning the beam 



physics is sufficient to address the closure issue. Specifically, we assume that the amplitude 
of the betatron oscillations in the frame moving with the beam is small compared to the 
characteristic size of the beam by the ratio 5, as we will discuss in more detail in Section 
III II when we describe the multiscale expansion. In this regime, we prove in the Appendix, 
starting from the collisionless Vlasov equation, that the pressure tensor P must be gyrotropic 
to lowest order: 

P =p±I+ 0|| -p±) bb + 0{5) (12) 

where b = B/5 is equal to e z in our geometry, p±_ is the pressure in the radial- longitudinal 
plane, and p\\ the pressure along the direction of the magnetic field. The divergence of the 
gyrotropic pressure tensor is easily calculated: 

V • P = Wp± + (p\\ - p±) b ■ Vb + b • V -p±)b+ On - P± ) V • bb + 0(5) (13) 

Now, since b • Vb = e z ■ Ve z = 0, in the radial-longitudinal plane the divergence of the 
pressure tensor takes the form of a pressure gradient to lowest order in 5: 

(V-P) ± = V P± + 0(5) (14) 

In the multiscale analysis presented in Section IHfl we will need to expand the beam fluid 
equations to second order in 5. Since in Eq. (TT0|) the pressure term is multiplied by S 2 , we 
only need an expression for the pressure tensor to zeroth order in 5. This is precisely what 
we have in Eq. (114jl . In the rotating frame, the appropriate fluid equations for the beam 
thus take the following form 

dn „ / \ 

z + nV-v = 15 
at 

^1 + v x ez = _ 5 2 ( V0 + ?L Vp \ (16) 
at \ n J 

V 2 = -n (17) 

where for the simplicity of the notation we have dropped the _L symbol for p since we only 
consider the two-dimensional radial-longitudinal plane and p\\ never enters in the analysis. 

Equations ffT5"|) - fTTT|) are the fluid equations we solve in the remainder of this article, using 
a multiscale analysis. Note that these equations are not fully closed in the sense that an 
equation for the evolution of the pressure p is missing. As we will show shortly, such an 
equation is not necessary. Indeed, when the pressure term can be expressed as an exact 
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gradient in the momentum equation, pressure effects do not enter in the final equation for 
the evolution of the beam density on the space charge time scale. 



III. MULTISCALE ANALYSIS OF THE FLUID EQUATIONS 
A. Multiscale expansion 

Given the difference in the time scale associated with the betatron oscillations and that 
associated with the expansion of the beam due to electrostatic and thermal forces, we perform 
a multiple time scale analysis of Equations (I15]) - (jl7p . Specifically, each quantity Q is assumed 
to vary according to the different time scales as follows: 

Q(r, t) = Q(r, t , * 2 , U, ...) = Q(r, t, 5% 5% . . .) (18) 

With this formal expansion, we have 

It is convenient for the rest of the calculation to separate the quantities Q into the sum of 
a rapidly oscillating part Q due to the betatron oscillations, and a slow monotonic evolution 
Q due space charge and thermal effects: 

Q(r, t , t 2 , ...) = Q(r, t Q , t 2 , . . .) + Q{r, t 2 , . . .) (20) 

Observe that, by construction, Q does not depend on to- We will show shortly that Q is 
periodic in t . 

In the frame moving with the beam, the first non-vanishing contribution to the velocity 
is due to the betatron oscillations. As mentioned previously, we assume in our analysis that 
the amplitude of these oscillations is of order 8a. Since the betatron time scale is u c , this 
means that the first non-vanishing contribution to the velocity is of order Sau c . In other 
words, the velocity comes in first order in 5 in our ordering and the appropriate expansion 
for the relevant physical quantities is the following: 

n = n + 5 (ni + ni) + 5 2 (n 2 + n 2 ) + 0(S 3 ) (21) 

p = p + O(5) (22) 

= 0o + O(5) (23) 

v = 5vi + 5 2 (v 2 + v 2 ) + 0(8 3 ) (24) 

10 



We introduce this expansion in the set of equations given by Eqs.( !T5|) - (fl7|) . and solve 
order by order in 5. We start with Poisson's equation, which we only need to zeroth order 

V 2 0o = -no (25) 

Turning to Eq. ( TT5|) . the mass conservation equation, we see that the first non-trivial equation 
arises in first order in 6 and is given by 

^± + V ■ (nbvx) = (26) 

This equation describes the evolution of the density on the fast time scale, due to the 
betatron oscillations. To next order, the mass conservation equation takes the form: 



dri2 dn 



V ■ [(ni + rii) vi + n (v 2 + v 2 )] = (27) 



dt dt 2 

Averaging this equation over the fast time scale, we get, because of the periodicity of the 
oscillating quantities in t : 

dn 



dt 2 

where we have introduced the notation 



V • (< rixVi > +n v 2 ) = (28) 



2tt 



< Q >= 77- / Qdt 
2tt Jo 

Eq.f l28p describes the evolution of the bunch density on the slow time scale. This is where 
the influence of the electrostatic effects and the beam temperature will enter, and this is 
therefore all we need from the mass conservation equation. The purpose of the remainder 
of the article is to solve this equation, first by deriving expressions for < niVi > and v 2 
in terms of zeroth order quantities, and by then solving the resulting equation numerically. 
The necessary expressions are obtained from the momentum equation. To first order in S, 
we have 

^+vixe 2 = (29) 
ato 

and by averaging the second order momentum equation over the fast time scale as was done 
before, we find 

a 2 

v 2 =< Vi ■ Vvi > xe z + V0 O x e z + — Vp x e z (30) 

n 

At this point, all the relevant equations have been derived. From Eqs. (l26|) . (129]) and 
fl30|) we can derive expressions for hi, vx and v 2 in terms of zeroth order quantities and 

11 



initial conditions, which we can then insert in Eq. (l28j) to obtain the desired equation for 
the evolution of the beam density on the space-charge time scale. This is done in the next 
subsection. 



B. Solving the equations 

We start the calculation with the description of the betatron motion to lowest order, i.e. 
vi, and its effect on the bunch density n\. Vi is given by Eq. (l29|) . which is just the equation 
of a particle immersed in a uniform magnetic field. It is easily solved: 

vi(r, t , t 2 , ■ ■ •) = uicosto + e z x uisint (31) 

where Ui(r, t 2 , . . .) is the initial betatron velocity, i.e. Vi at time to = 0, 2ir, . . .. Now that vi 
is known, Eq. (T2"6]) is most easily solved by introducing the displacement vector ^ 1 (r, to, h, ■ • •) 
defined by 

vi = ^; li(r,t = 0,t 2 ,...) = (32) 



Substituting for vi with £± in Eq. (l26|) . we find 

d_ 

which is easily integrated as 



nx + V-lri&A =0 (33) 



Tii = Ni (r) - | x ■ Vn - n V ■ | x (34) 



where N\ is the initial beam density. Eq. fl32|) is integrated as easily using Eq. (|3ip . and we 
have 

£i = uisinto + ui x e z (costo — 1) = (vi — Ui) x e z (35) 

These results to first order in 5 can now be used in Eq. (1301) and Eq. (j28p . After a slightly 
lengthy calculation, we find 

< niVi > +n < vi • Vvi > xe 2 = x (fiQu\e z ) (36) 

so that the equation for the evolution of the beam density on the space-charge time scale is 
given by 

^ + V • (n o V0 o x ej + a 2 V • (Vp x e 2 ) = (37) 
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From the identity Vpo x e z = V x (po^z) we see that the pressure term vanishes in Eq.([3 
and we finally have 

^ + V0 O x e z • Vnb = (38) 

012 

Eq. (l38p . combined with Poisson's equation 

V 2 0o = -no (39) 

forms the desired closed set of coupled equations, describing the evolution of the beam 
density under the effect of space-charge forces. It is interesting to note that temperature 
effects do not enter in this equation when the pressure tensor is gyrotropic, as it is with 
the ordering we chose for the amplitude of the betatron oscillations. This is the reason why 
we do not need to add an equation for the evolution of the pressure to the fluid equations 

(HSD-dHD- 

Eq. fl38|) can be easily interpreted as the advection of the density profile in the velocity 
field E x B/5 2 , the so-called E x B velocity (there is no - sign in front of V0o because 
brce is in the opposite direction in the moving frame). This is in agreement 



the Lorentz 

with Eq.(l) [20|, and can be explained as follows. In the absence of accelerating gaps, the 
effect of the betatron oscillations on the density profile averages out to zero on the slow time 
scale. Thus, the density profile simply follows the slow, averaged motion of the ions. In a 
homogeneous magnetic field, in the presence of an electric field created by the charges in 
the bunch themselves, this slow motion is just the E x B motion. 

With this intuitive interpretation, it is easy to predict the early stages of the evolution of 
the beam for the density profiles discussed in the Introduction. The two situations of interest 
are illustrated in Figures |3] and HI which show the V</>o x e z velocity field due to a cylindrically 
symmetric density distribution and an elongated density distribution respectively. If we 
start, as in Figure [31 with a density profile which is cylindrically symmetric, i.e. such that 
the beam density is only a function of the distance from the center of the bunch (n = no(r)), 
then the initial electrostatic potential due to the charges in the beam will also be cylindrically 
symmetric: 0o = 0o( r )- This implies that the initial electric field will be along the e r 
direction, E = E r e r , and the convective velocity V0o x e z will be along the e# direction. We 
thus expect the density profile to rotate around the center of the bunch. Since this density 
profile is cylindrically symmetric, a rotation of the profile around the center will leave the 
profile invariant. In other words, the space charge forces in a cylindrically symmetric bunch 

13 



in a purely drifting region are such that the bunch pro per ties are kept constant. This result 



21] with a different method and with 



agrees with a similar result obtained by Kleeven et al. 
a very different set of assumptions. 

If, however, the initial density distribution is elongated in the direction of propagation of 
the beam as in Figure HJ the cylindrical symmetry will be broken, and the E x B velocity 
field will distort the beam. Since the electric field is strongest at the extremities of the 
beam, the distortion will be strongest in these places. This initial situation explains the 
formation of the spiral arms in the early stages of the beam evolution. In addition, since the 
electrostatic potential far away from the beam still retains cylindrical symmetry (far away, 
the beam is seen as a point charge), we expect the E x B convection to act in such a way 
as to make the beam conform with this cylindrical symmetry. This is the reason for the 
formation of the often mentioned galaxy-like shape, and after more turns, the formation of 
a compact, stable cylindrically symmetric core. In the next section, we solve the system of 
coupled equations Eq.( l38|) and Eq.(|39l) numerically and confirm these predictions. 



IV. NUMERICAL RESULTS 



In order to solve Eq. (1381) and Eq. (l39l) numerically, we first rewrite them as follows 

dn 2 f dn dn 



<«, 

In Equation ()40|) . t represents the numerical time steps. Since they are defined in terms of lo c 
in the simulations, i.e. in terms of the fast time scale, 5 2 appears in the advection term. E x 
is the electric field in the radial direction, satisfying E x = —d<f)o/dx and E y is the electric 
field in the longitudinal direction, satisfying E y = —d(j)o/dy. We assume open boundary 
conditions for the electric field, so that E is simply due to the charges in the beam. 

We solve Eq. (I40j) and Eq. (1411) by discretizing the problem on a regular Cartesian grid. The 
time integration of Eq. (l40p is done with a second order backward difference scheme, except 
for the first time step which is calculated with a second order trapezoidal scheme. At each 
time step, we compute the partial derivatives of the density at each grid point with second 
order centered differences, assuming that the beam density is zero outside the boundary of 
the computational domain, which is always chosen to be larger than the maximum beam 

14 



-2 



-1 1 

Radial direction 



2 



FIG. 3: Velocity field V0o X e- (magenta arrows) for a cylindrically symmetric density distribution. 



To compare with figures in Ref. 
y direction. 
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and Ref. 



161 ] . the beam's transport direction is along the negative 



size. We also compute E x and E y at each time step on this grid by numerically solving the 
two Poisson's equations in Eq. (l4ip . Specifically, we calculate E x and E y at interior points 
by inverting the standard finite difference second order Laplacian matrix, and use the two- 
dimensional free space Green's function for Poisson's equation to calculate the electric field 
at the computational boundary: 

») = s J J dx 4y { —^- I — w (42) 

EM, ») = Y ,]\ dx d * fa-^ + te-yT (43) 



In Eq. (l42p and Eq. (1431) . the subscript b stands for points on the boundary of the computa- 
tional domain. We use Simpson's rule to perform the numerical integration. Combining all 
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FIG. 4: Velocity field V^>o x e z (magenta arrows) for an elongated density distribution. To compare 



with figures in Ref . [3| and Ref . 



161 ]. the beam's transport direction is along the negative y direction. 



these elements, we have a numerical scheme to integrate Eq. (140!) and Eq. (|4"TI) that is second 
order in both time and space. 

For the simulations we present in this article, we choose the same beam parameters as 

n 

the ones used in [3[ for the study of a 1mA, 3 MeV coasting beam in PSI injector II. With 
these parameters, we find 5 2 ~ 0.8. This value is at the upper limit of the validity of the 
multiscale analysis used in this work, but will allow us to compare our results with a large 
300!^ of PIC simulations performed to investigate space charge effects in PSI injector II 
9]. In agreement with sj], we start with an initial charge distribution given by 



n(x, y) = iVoexp 



x 



2o* 



(44) 



1(3 



A. Stability of a round beam 

We first demonstrate numerically that a round beam is completely stationary under the 
effect of space charge forces, as we discussed in the previous section. We choose 2a x = 2cr y = 
2.52 mm and run our simulation until the beam has completed 40 turns. The results are 
shown in Figure [5j It is clear in this figure that the charge distribution is not affected by 
the advection in the V0o x e z velocity field, and the charge distribution after 40 turns is 
identical to the initial charge distribution, as expected. In fact, the results shown in Figure 
[5] can be seen as a proof of the accuracy of our numerical method. We can now confidently 
turn to the study of elongated beams and beam spiraling. 



o 



Radial direction 



I 01 



-5' 

-5 



o 



Radial direction 



-5" 

-5 



o 



Radial direction 



n 



I 



0.8 
0.6 
0.4 
0.2 




-5' 

-5 



o 





Radial direction 



a 01 



-5 1 

-5 



o 



Radial direction 



5 1 

-5 



o 



Radial direction 



I 



0.8 
0.6 
0.4 
0.2 





FIG. 5: Top view of a round coasting beam in the local frame moving with the beam center. Up: 
turn 0, 5, 10. Down: turn 20, 30, 40. To compare with figures in Ref.[s]] and Ref. [if]], the beam's 
transport direction is along the negative y direction. S 2 = 0.8 and 2a x = 2a y = 2.52 mm. 
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B. Elongated beams and beam spiraling 



For the study of beam spiraling, we consider an elongated beam with 2a x = 2.52 mm and 
2a y = 13.4 mm (15° phase width) as in jjj. The numerical results are shown in Figure 
(early times) and Figure [7J (later times), where, as in Ref. [sj] and Ref. [if]], we have truncated 
the beams at 10% of the maximum charge density. 



By comparing Figure [6] with the early stages of the beam evolution shown in Ref. 16|, we 
observe that there is a strong similarity between our results and those obtained in simulations 
with the PICS code. We obtain a slightly stronger deformation than seen with PICS in 
Ref. [3]. This small discrepancy can be attributed to two factors. First, the magnetic 
geometry is different in the two simulations: the magnetic field in our simulation is uniform 
and in the z direction, while Adam used an analytic approximation of the magnetic field in 
the PSI Injector II cyclotron. Second, the shape of the initial charge distribution, which we 
chose to be identical to that in jjj, is slightly different. The good agreement between our 
results and the PICS simulations is not surprising since PICS is a two-dimensional PIC code 
which relies on the separation between the betatron and space charge time scales to reduce 
the computational complexity of the PIC simulations. The agreement also suggests that the 
details of the magnetic geometry do not play a strong role in the radial-longitudinal vortex 
motion due to space charges. 

Comparing Figures [6] and [7] with PICN results [3] and with OPAL-CYCL results 3], we 
conclude that there is qualitative agreement between all descriptions, with the formation 
of a compact stable core after 40 turns. Furthermore, as in jjj, we see that a tenuous low 
density halo persists, even after 40 turns. However, there remains small differences between 
[31 and our study. First, the core that is formed after 40 turns is round in our simulations, 

n 

unlike the slightly elongated core observed in [3J]. Second, the beam spiraling effect appears 
to be stronger in our results. As before, part of the discrepancy can be explained by the 
difference in the magnetic geometry assumed in the different simulations. It is also a direct 
consequence of the two-dimensional description that we have adopted in our study, which 
leads to a stronger self-electric field, and thus a larger V0o x e z velocity field. This effect has 
already been observed in the difference between the PICS and PICN codes Obtaining 
fully quantitative agreement with OPAL-CYCL on this problem would thus require a three- 
dimensional fluid treatment, and the inclusion of the details of the PSI injector II magnetic 
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FIG. 6: Top view of an elongated coasting beam in the local frame moving with the beam center. 
Up: turn 0, 1, 2. Down: turn 3, 4, 5. To compare with figures in Ref.Q| and Ref.ljjJ], the beam's 
transport direction is along the negative y direction. 5 2 = 0.8, 2a x = 2.52 mm, and 2a y = 13.4 mm 

geometry. 



V. CONCLUSION 



In the vast majority of high-intensity cyclotrons and rings, the betatron time scale remains 
much shorter than the time scale associated with space-charge and temperature effects. This 
separation of scales can be used advantageously in analytical and numerical studies of space- 
charge effects on beam transport. In this article we performed a multiscale analysis relying 
on this scale separation to study the radial-longitudinal vortex motion of an intense charged 
particle beam in a regime in which the pressure tensor is gyrotropic. As a result of this 
analysis, we have been able to show that the evolution of the charge distribution on the 
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FIG. 7: Top view of an elongated coasting beam in the local frame moving with the beam center. 
Up: turn 0, 5, 10. Down: turn 20, 30, 40. To compare with figures in Ref.^] and Ref.flf]], 
the beam's transport direction is along the negative y direction. S 2 = 0.8, 1a x = 2.52 mm, and 
2a y = 13.4 mm 

space-charge time scale is determined by the nonlinear advection equation of the charge 
distribution in the E x B velocity field, where E is the self electric field due to space charges 
and B is the externally applied magnetic field. The stability of intense round beams is easily 
understood in this light. A beam with a charge distribution which is only a function of the 
beam radius has an associated E x B velocity field which is always perpendicular to the 
beam radius in the radial-longitudinal plane. Such a velocity field cannot modify the charge 
distribution. Beam spiraling in elongated beams is explained just as easily: in such beams 
the cylindrical symmetry of the charge distribution is broken, and the E x B velocity field 
therefore modifies the charge distribution until the cylindrical symmetry is restored. The 
spiral shape is a natural consequence of the fact that E is strongest at the extremities of the 
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beam. 

Aside from the shape of the beam, the parameter that determines the strength of the 
spiraling effect is 5 2 = ujp/u 2 . In cyclotrons and rings with very intense beams and moderate 
magnetic fields, this parameter is relatively large (5 < 1) and beam spiraling can play a 
positive role since a dense round stationary core can be formed after a few turns and the 
spiral halo can be stripped at low energies. In accelerators with less intense beams or higher 
fields, however, such that 5 <C 1, beam spiraling can have a negative effect since it takes 
more turns to develop, and high energy beam halos can be formed in this process. 

The good agreement between our numerical simulations and PIC simulations demon- 
strates that the potential coupling between the radial-longitudinal motion and the vertical 
motion is not a key physics elements that need to be included to understand space charge 
effects in the radial-longitudinal plane. This suggests that the simplified fluid equations 
Eqs JToTrTTt together with an equation for the evolution of the fluid pressure p, might be very 
well suited for numerical studies of the novel and interesting regime in which 5 2 > 1. In 
this regime, of interest in some new machines such as the UMER ring 19], the analytic 
multiscale expansion used in this work breaks down, but the 2-D fluid equations presented 
here could still represent a significant reduction in complexity in numerical simulations as 
compared to PIC simulations. This topic is currently being investigated, with progress to 
be reported at a later date. 



Appendix 

In this appendix we prove, starting from the collisionless Vlasov equation, that the pres- 
sure tensor P is gyrotropic to lowest order in 5 when the amplitude of the betatron oscilla- 
tions is small by 5 compared to the size of the beam. 

The collisionless Vlasov equation determining the evolution of the beam distribution 
function f(r,v,t) is 

df e 

+ v • V/ + - (E + v x B) • V v / = (A.l) 
at m 

Let us compare the relative strength of each term in this equation relative to the magnetic 
term, in the regime where v ~ 5au c : 

E ena/e uj 2 mv • V/ 5oo c mdf jdt 5u c 

v x B 5au 2 m/e 5u 2 ev x B • V v f u c ev x B • V v / oj c 
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We see that all the terms are small by 5 compared to the magnetic term, so that to lowest 
order in 5, Eq. flA.lj) is: 

vxB- V v f = + 0(5) (A.2) 

Equation (1A.2|) implies that to lowest order in 5 the distribution function 
f(x, y, z, v x , v y , v z , t) has the following particular dependence on v x and v y : 

f(x, y, z, v x , v y , v z , t) = f(x, y, z, \Jv 2 x + v%, v x , t) + 0(5) (A.3) 

Thus, to lowest order / is even in v x and v y , which implies that to the same order the 
off-diagonal components of the pressure tensor vanish: 

P = Pil + (p\\ - P±) bb + 0(5) (A.4) 
with p ± = p xx = p yy and p\\ = p zz . 
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